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Abstract 

Equilibrium equations in the form of complementarity conditions often appear as 
constraints in optimization problems. Problems of this type are commonly referred 
to as mathematical programs with complementarity constraints (MPCCs). A popu- 
lar method for solving MPCCs is the penalty interior-point algorithm (PIPA). This 
paper presents a small example for which PIPA converges to a nonstationary point, 
providing a counterexample to the established theory. The reasons for this adverse 
behavior are discussed. 
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1 Introduction 

Equilibrium equations in the form of complementarity conditions often appear as con- 
straints in optimization problems. Problems of this type are commonly referred to as 
mathematical programs with complementarity constraints (MPCCs), and arise in many 
engineering and economic applications; see the survey [7j and the monographs ^HJ QB] for 
further references. The growing collections of test problems [T4"l U] indicate that this an 
important area. MPCCs can be expressed in general as 

minimize f(x,y,w,z) 

subject to x E X . N 

7-1/ X n 1-1 

b [x, y, w, z) = 
< y JL w > 0, 

where X C IR n is a polyhedral set, / and F are twice continuously differentiable functions, 
and y,w G IR m . The complementarity constraint y _L w means that either a component 
of y is zero, or the corresponding component of w is zero, which implies that y T w = 0. 

Many approaches exist for solving MPCCs. These include branch-and-bound methods 
|2j, implicit nonsmooth approaches [T%] . piecewise SQP methods [IB], perturbation and 
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penalization approaches analyzed in [22], and the penalty interior-point algorithm 
(PIPA) [Oil, whose convergence properties are the subject of this note. 

Recently, there has also been renewed interest in solving (jl.ljl with standard nonlinear 
programming (NLP) techniques by replacing the complementarity constraint by 

y T w < 0, y > 0, w > 0. 

This approach had been regarded as numerically unsafe because fll.lj) violates a constraint 
qualification at any feasible point |21j . However, numerical [B] and theoretical evidence 
[TJ IS] suggests that sequential quadratic programming (SQP) solvers are able to solve 
large classes of MPCCs. A similar analysis is being extended to interior-point methods 
(IPMs) 0II3l2n]- 

This note is organized as follows. The next section briefly describes PIPA and presents 
its convergence result. Section |3] presents a small counterexample that shows that PIPA 
can fail to converge for certain admissible choices of parameters. Finally, Section 0] dis- 
cusses the example and suggests a remedy that avoids this problem. 

2 The Penalty Interior-Point Algorithm (PIPA) 

PIPA PH Chapter 6.1] solves MPCCs of the form ([Ol. It has been studied in [TT1 ITH] 
and has been applied to a number of real applications in finance [13 , target classification 
[T7] . and electricity markets [T2] . 

PIPA generates a sequence of iterates (x k , y k , w k , z k ) that are strictly interior in y k > 0, 
w k > 0. The iterates are computed by solving a quadratic direction-finding problem. In 
that sense, PIPA combines aspects of interior-point and SQP methods. 

We use d to denote the step or displacement computed by PIPA. Subscripts such as d x 
refer the part of d corresponding to the x- variables; superscripts are used to denote iterates 
or evaluation of functions at a particular point, for example, V/ fc = V f(x k ,y k ,w k , z k ). 
Diagonal matrices are denoted by W = diag(u>) and Y = diag(?/) and e = (1, . . . , 1) T . 

The algorithm solves the following direction-finding problem at every iteration: 

minimize V ' f kJ d + \d T x Q k d x 
subject to x k + d x E X 

F k + VF kT d = ( 22 ) 

Y k d w + W k d y = -Y k w k + (J^^e 

Kill < c(j|F fe || +y kT w k ^j , 

where c > and a G (0, 1) are parameters and Q k approximates second-order information 
in / and F. Note that (|2.2j) has a quadratic trust-region-type constraint. If this constraint 
is replaced by the norm, then ()2.2j) becomes a quadratic program (QP). Note also that 

the linearization of the complementarity constraint Yw = is relaxed by a v ™ e, giving 
PIPA an interior-point flavor. Note that PIPA is not strictly an interior-point algorithm, 
as it remains interior only with respect to the complementary variables y, w but not with 
respect to the other constraints. We will show that the final trust-region type constraint 
is the source of the problems affecting PIPA. 
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After determining a search direction, PIPA performs a backtracking line-search to 
ensure that y k , w k remain positive and close to the central path. This is achieved by 
finding the root r of 

v kT w k ( d kT d k \ 

9k (r) = (1 - p)o V -^- + r min « - ■ (2-3) 

m \ i<i<m y ' ' ml 

Global convergence is enforced by further reducing the step length to yield sufficient 
reduction in a quadratic penalty function 

P a (x, y, w, z) = f(x,y,w,z) + a(\\F(x,y,w,z)\\ 2 + y T w) , (2.4) 

where the penalty parameter a > is chosen to be the smallest a such that the step 
ensures sufficient model decrease, 

Vf T d - a p {l - a)y kT w k < -a p (l - a)y kJ ' w k < -y kT w k . (2.5) 

Note that the penalty function P a mixes a quadratic penalty for the nonlinear equation 
F with an £± penalty for complementarity (since y k ,w k > 0, it follows that \\y k w h \\i = 
yk w ky Further details of the algorithm can be found in [TS|. PIPA can be summarized 
as follows. 

Penalty Interior-Point Algorithm (PIPA) 

1. Initialization: 

Choose parameters c > 0, a, 7, p G (0, 1), 

choose a starting point y°, w°, z°) such that y°, w° > suitably centered, 
choose a penalty parameter a > and set k = 0. 
REPEAT 

2. Direction finding problem: Solve problem (|2.2|) for a trial step d = d k . 

3. Step size determination: Find a step size r = to: 

3.1 Ensure centrality and positivity of (y,w) by finding the root of </fc(r) in ()2.3|) 
or setting r = 1 — e if this root does not exist or is greater than 1. 

3.2 Ensure sufficient reduction in the quadratic penalty function P a (x, y, w, z) 
in (12 .4j) by performing an Armijo-search on P a . 

Let 77 be the step size determined in 3.1 and 3.2. 

4. Update: (x k+1 , y k+1 , w k+ \ z k+l ) = (x k ,y k ,w k , z k ) + r k (d k ,d k , d k w ,d k ), k = k + l. 
UNTIL < e 

This presentation of PIPA is less sophisticated than ^B], which allows for instance a 
to vary with the iteration. Convergence of PIPA is established in [TH] under the following 
two assumptions. 

[SC] Strict complementarity of the solution, namely y* + w* > 0. 

[NS] Nonsingularity of the matrix 

" \/ y F* \7 W F* V Z F* ' 
W* Y* 
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A sufficient condition of [NS] is that the Jacobian of F satisfies the mixed P® property. 
The following theorem summarizes the convergence results in |16j . 



Theorem 2.1 [16, Theorem 6.1.17] Suppose that the Hessian matrices W k are 
bounded and that a > 0. If the penalty parameter a is bounded, then every limit 
point of (x k , y k , w k , z k ) that satisfies [SC] and [NS] is a stationary point of ([Lip . 



In the next section, we present a small example for which PIPA converges to a non- 
stationary point, contradicting Theorem 12.11 

3 A Counterexample 

In this section, we show that PIPA may fail to converge to a stationary point. This result 
contradicts the convergence results of [16j. This section examines the behavior of PIPA 
applied to the following small example: 

minimize x + w 

subject to — 1 < x < 1 , , 

-l+x+y=0 \ J 

< y JL w > 0. 

The solution of (j3.1j) is (x, y, w)* = (—1,2,0), which satisfies the assumptions [SC] and 
[NS], since y* + w* = 2 > and 



VyF* 


V F* ' 

v w 




' 1 


" 


w* 


Y* 







2 



is nonsingular. Next, we show that this example generates a limit point that is not 
stationary. Specifically, we perform a standard iteration of the algorithm and bound the 
limit away from stationarity. The following lemma summarizes our result. 



Lemma 3.1 PIPA applied to the Example ()3.1|) and started at (x°,y°,w°) = 
(0, 1, 0.02), with c = 1, a — 0.1, 7 = 0.01, p = 0.9 and a = 2 generates a sequence of 
iterates satisfying 





1 < y k < y k+l 


(3.2) 




1 < l w k < 2 
- 2 ~ 100 


(3.3) 


x k+1 > x k 


— a/ y k w k > —1 


(3.4) 


y 


" V +1 < -y k w k 
~ 2 y 


(3.5) 




5 

- < Tfc < 1. 

9 - 


(3.6) 
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Proof. The starting point satisfies the linear equation, which implies that they are 
satisfied for all subsequent iterations. The theorem follows by induction. Clearly, the 
assertions hold for the starting point. Now assume that ()3.2)1 to ()3.6j) hold for k — 1, and 
show that they also hold for k. 

The direction-finding problem for this example can be simplified as follows: 



minimize d x + d w 
subject to — 1 < x k + d x < 1 
d x + d y = 

w k d y + y k d w = (a — l)y k w k 
-s/y k w k <d x < ^/y k w k , 



(3.7) 



where Q k = for simplicity. However, the example remains valid for certain positive 
semi-definite bounded Q k . 

First we show that the optimal basis B to the LP (|3.7|) is given by 



" 1 





l " 


-1 


" 


1 


1 


w k 
















y k 







1 


-1 



From this, multipliers can be computed as 



B~ 








y 



o 4 

y 

i -i 4 

yk 




/ -4 \ 



Optimality follows, since A3 > for the only inequality constraint. The step d can also 
be computed as 

\ 

d k 



-Tuk 



B- J b 



/ 

V 



\/y k w k 



10 w 



where b k = (0, — j^y k w k , —\Jy k w k ) T is the right-hand side of the active constraints. 

Note that for any positive semi-definite Q k , this solution remains unchanged as long 
as A3 > 0. The contribution of the Hessian to the multipliers can be estimated as 



ill 

\ k + B~ 1 Q k d k > 1 + — - llB-^lllld* 1 



which remains positive for ||<5 fe || or \\d k \\ sufficiently small. Thus, the same conclusions 
apply even if a positive definite Q k is used. 
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Next, we show that (|3.4|) follows by applying the step d k , using r k < 1 and induction: 



x 



k+1 



x k + r k d 



kjk 



x 



> x k 



> x k 



y/y~w 



fc„„& 



-y/2 
10 



z=o 
/ 



y°w° 



1 



(>/2 



V 



1-4? 



v/2 



(3.f 



/ 



It is easy to show that the expression in the last line is greater than — 1, which implies 
that the "choice" of active set in the LP ()3.7j) was correct, and the lower bound on x is 
never active during the iteration. Next we show that (|3.2j) follows by applying the step 
in the y component, 

y k+l = yk + T k^ y k w k >y k >\. 

Next consider (|3.6|) . and observe that r k < 1 follows trivially. To obtain the lower bound, 
observe that the root of ^fe(r) in (|2.3J) is given by the following expression, where the 
superscripts k have been omitted for the sake of simplicity 

,2 



y 



(9y + lOy/yw) ^fyw 



1 yw 
10 d w d y 

1 y/V 1 1 y/V 

V ^ 9 V 



E 

1=0 



10 



9 



From the fact that w < jjjg, y > 1, induction, and > 0, the step size can be 



bounded below by its leading terms: 



729 yjy 



T > 



1 1 



> 10_L 10 



> 



9 ^2 
5 

9 ' 



81 



10 1 100 Vw 
81 ^ 729 ^ 
20 
6561 



1000 



6561 1 



+ 



It remains to show that this lower bound on r also passes the second line-search criterion. 
The penalty update rule in (6.1.24)] finds the smallest integer p > 1 such that 



Vf kT d - a p (l - a)y k 'w k < -a p (l - a)y 



k T k k T k 

s w k < -y k w k 



holds, where we have used the fact that ||F fe || = for all k. Since Vf k d < and 
a(l — a) > 1, this is always satisfied for p — 1, and the penalty parameter is never 
increased. Finally, this step size also satisfies the conditions for the Armijo search in Step 
3.2. The actual reduction clearly satisfies 



Since y 



k+l w k+l 



Pa +1 ~ Pa = rd x + rd w + a (y 
< \y k w k , it follows that 



k+1 w k+1 



- y k w k ) . 



P k+l - P£ < rd x + rd w - a-y*w 



k„..k 
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which implies the sufficient reduction condition, 

Pa +1 ~Pl< ird x + ^Td w -a{l- a) iy k w k , 

since 7 < 1 and (1 — a) r y < ~. 

Next consider ()3.3|) . Since r > |, it follows that 



w 



k+1 



9 k ™ r~k k\ / 1 k 5w r~k I 

— w rVv w 1 — ~ w rVy w ■ 

10 y k / 2 9 y k 



w + t \—- 

Thus, it is possible to bound the new complementarity violation after the step by 



(3.9) 



y 



y k w k 



l„,,fc 5 w 
9 y k 



< ±y w , 



9 



w 



ln,K„k _ 5„,k^ y k w k + l w k^ y k w k _ ji^fc) 



k\2 



obtaining ()3.5|) . This concludes the proof. 



□ 



By taking the limit (A; — > 00) in ()3.8|) . it follows, that all control iterates satisfy the 
bound 

x k+1 > — « -0.4828. (3.10) 

~ 10(V2-1) 

Similarly, an upper bound can be derived for the state iterates as 



,k+l 



1=0 



V2 
10 



1 \V2 



1 " -4 



which shows that 



,k+l 



< 1 



10 



1) 



1.4828. 



(3.11) 



These bounds can also be verified numerically. Example ()3.1|) and PIPA have been im- 
plemented in AMPL PHI by using the "looping extension", which allows the convenient 
implementation of algorithms in AMPL. The sequence that is generated is given in Tabled 
and confirms the bounds given above. The columns headed "ared" and "pred" in Tabled 
refer to the actual and predicted reduction. Their implications are discussed below. 
Thus, it follows that all iterates remain in a compact set, namely, 




1 



x° > x k > -0.4828 
y° < y k < 1.4828 



m = w° > w k > 0, 



and the sequence (x, y, w) k has a limit point (x, y, w)°°. From ()3.9|) it follows that 



w k+1 < w k 



w 
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k 


h 

X 


yk 


k 

w 


ared 


pred 


1 





1 


0.02 






2 


-0.096022613 


1.0960226 


0.0058578644 


-0.198 


-0.137 


3 


-0.17606958 


1.1760696 


0.00016323495 


-0.0974 


-0.0982 


4 


-0.18991126 


1.1899113 


1.4549224E-05 


-0.0143 


-0.0143 


5 


-0.1940679 


1.1940679 


1.4171928E-06 


-0.00421 


-0.0042 


6 


-0.19536745 


1.1953675 


1.4145236E-07 


-0.00131 


-0.0013 


7 


-0.19577825 


1.1957782 


1.4223933E-08 


-0.000412 


-0.000411 


8 


-0.19590853 


1.1959085 


1.433645E-09 


-0.00013 


-0.00013 


9 


-0.1959499 


1.1959499 


1.4460519E-10 


-4.14e-05 


-4.14e-05 


10 


-0.19596304 


1.195963 


1.4586827E-11 


-1.32e-05 


-1.31e-05 



Table 1: Iterates, actual and predicted reduction of AMPL implementation of PIPA 



so that w°° = and y°° w 00 = 0. Therefore, PIPA converges to a limit point 
such that 

> x°° > -0.4828 

1 < y°° < 1.4828 

w°° = 0. 

This limit point satisfies complementarity and the linear constraints for problem (|3.1|) 
(since (x,y,w)° satisfies the linear constraints). Moreover, this limit point also satisfies 
the assumptions of Theorem 12 .11 Assumption [SC] is satisfied because y^ + w 00 > 1 > 0, 
and Assumption [NS] holds because 



" M 


I 




"10" 


w°° 






y°° 



is nonsingular for all 1 < y°° < 1.4828. However, the limit point is clearly not a stationary 
point, thereby contradicting Theorem 12.11 

4 Discussion and Remedy 

In the preceding section, a small example has been presented for which PIPA converges to 
a feasible but nonstationary limit point. This example contradicts ^Sl Theorem 6.1.17], 
which establishes the convergence of PIPA. The example might be considered contrived 
in the sense that we choose a particular starting point and parameter values. However, 
for any different starting point it is possible to construct parameters such that PIPA fails. 
The reason for this apparent failure of PIPA is the trust-region type constraint 

Kill < c(\\F k \\+y kT w k ) =: A k (4.12) 

in the direction- finding problem ()2.2j) . The trust-region radius converges to zero as 
the iterates approach feasibility, thereby limiting the progress toward optimality in the 
controls x. This adverse behavior can be expected to occur whenever the iterates approach 
feasibility "faster" than optimality. 
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Linking the trust-region radius to the feasibility of the problem is counterintuitive. 
Normally, A& is controlled by the algorithm, rather than by the iterates. In particular, 
it should reflect how well the model problem ()2.2j) approximates the original problem, 
measured, for instance, by the agreement between actual reduction, 

ared = P k a - P k+ \ 

and predicted reduction, 

pred = Vf kT d-a p (l-a)y kT w k . 

Table ^ shows that in this example, both the actual and the predicted reductions show 
near perfect agreement toward the end. Thus, one would expect Afc to be increased rather 
than decreased. The adverse situation in which (|4.12j) limits progress toward optimality 
can be easily detected. The Lagrange multiplier of (J4.12j) indicates whether progress can 
be made by relaxing this constraint. In the present example, this multiplier converges to 
1. 

One possible way to remedy the adverse behavior of PIPA is to control the trust- 
region radius by using standard trust-region techniques. The direction-finding problem 
then becomes 

minimize Vf kT d + \<$ L Q k d x 
subject to x k + d x G X 

F k + VF kT d = (4.13) 

Y k d w + W k d y = -Y k w k + a y -^e 
\\d x \\ < A fc , 

where the trust-region is now controlled by the algorithm. A typical updating scheme for 
Afc is as follows. First compute the ratio of actual over predicted reduction 

ared 

Pk ■= 7, 

pred 

and then update Afc according to 

{[70 A fc , 7! A k ] if pfc < rji 
[71 A fc , Afc] if p k G [771,%) 
[Afc, 72 Afc] if pfc > r] 2 

in Step 4, where < 70 < 71 < 1 < 72, and < 771 < 772 < 1. Typical values for these 
constants are 70 = 7i = 0.5, 72 = 2, rji = 0.25, and r] 2 = 0.75. 

This results in a mixed trust-region/line-search algorithm. Convergence can be estab- 
lished under the assumption that the Jacobian of F, 

' V y F V W F V Z F ' 
W Y 

remains nonsingular and has a bounded inverse along the lines of the coupled trust-region 
approach of Dennis et al. [3] , since the QP subproblem can be regarded as a problem in 
d x only after eliminating (d y , d w , d z ) using the linearization of F. 
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